#!/usr/bin/env python

""" Try to obtain the epsg code from a shapefile (that includes a .prj)"""

import os
import sys
import glob
import string

sys.path.append(os.path.dirname(__file__))

try:
    import termout
except ImportError:
    sys.exit("Requires script from http://stackoverflow.com/questions/1396820/ to be named 'termout'")
    
from osgeo import osr


def get_spatial_reference(shapefile):
    srs = osr.SpatialReference()
    shape_name = os.path.splitext(shapefile)[0]
    shpfile = os.path.basename(shape_name)
    prj_string = open(shape_name +'.prj','r').read()
    ret = srs.SetFromUserInput(prj_string)
    proj4 = srs.ExportToProj4()
    if not proj4:
        print "Failed to convert prj of %s, trying again using 'ImportFromESRI' method..." % shpfile
        prj_file = open(shape_name +'.prj','r')
        prj_lines = prj_file.readlines()
        prj_file.close()
        for i in range(len(prj_lines)):
            prj_lines[i] = string.rstrip( prj_lines[i] )
        srs = osr.SpatialReference()
        srs.ImportFromESRI(prj_lines)
        proj4 = srs.ExportToProj4()
        if not proj4:
            print "Failed to convert prj of %s, giving up..." % shpfile
        else:
            print "Second try assuming ESRI wkt worked for %s!" % shpfile           

        srs.from_esri = True
    else:
        srs.from_esri = False
    srs.AutoIdentifyEPSG()
    return srs

def to_proj4(srs):
    return srs.ExportToProj4()

def to_epsg(srs):
    if srs.IsGeographic():
        return srs.GetAuthorityCode('GEOGCS')
    else:
        return srs.GetAuthorityCode('PROJCS')

def shortname(srs):
    if srs.IsGeographic():
        return srs.GetAttrValue('GEOGCS')
    else:
        return srs.GetAttrValue('PROJCS')

def clip(str_, len_):
    if not str_:
        return str_
    if len(str_) > len_:
        return str_[:len_] + '...'
    return str_

def get_results(shape):
    srs = get_spatial_reference(shape)
    return str(to_epsg(srs)), str(srs.from_esri), str(clip(shortname(srs),20)), str(clip(os.path.basename(shape),15))


if __name__ == '__main__':
    
    print '\n' 
    
    table = []
    if not len(sys.argv) > 1:
        sys.exit('Please supply path to shapefile or a directory containing shapefiles')

    shape_or_dir = sys.argv[1]
    
    row = ('epsg','from esri','name','shape',)
    table.append(row)
    
    if os.path.isfile(shape_or_dir):
        table.append(get_results(shape_or_dir))
    
    elif os.path.isdir(shape_or_dir):
        shps = glob.glob(os.path.join(shape_or_dir,'*.shp'))
        if not shps:
            shps = glob.glob(os.path.join(shape_or_dir,'*.SHP'))
        for shape in shps:
            table.append(get_results(shape))

    print '\nSummary:' 
    
    termout.colprint(table)
